function [mm,opt_hp_N,cdf_N,cdf_O,avg_diff_N,avg_diff_O,policy_change_table] = DO_Appendix_D2(param_D2,psi)
    %% ==================== Set Parameters and Grids ========================%

    %----- Specify Parameters -----%
    nu = param_D2(1);
    beta = param_D2(2);            % Discount factor (factors in exit probability)

    alpha = param_D2(3);           % Production function parameter
    d = param_D2(4);               % Human capital depreciation rate
    ell = param_D2(5);             % Matching function parameter
    lambda = param_D2(6);          % Probability of search OTJ
    chi = param_D2(7);             % Investment cost parameter
    b=param_D2(8);                 % Unemployment benefit

    muA =param_D2(9);              % Staffing agency fee
    c =param_D2(10);               % Cost of posting a vacancy
    delta =param_D2(11);           % Separation rate
    h_sigma =param_D2(12);         % Standard deviation of initial human capital distribution

    pr_phi(1) = param_D2(13);      % Probability of entering the model with phi=1 (low fixed productivity type)
    pr_phi(2) = param_D2(14);      % Probability of entering the model with phi=2 (high fixed productivity type)

    z(1) = param_D2(15);           % Productivity parameter f(phi,h) = z(phi)h^(alpha)
    z(2) = param_D2(16);

    gamma=2;
    vfi_toler = 0.0001;
    dist_toler = 0.0000001;
             
    h_mu = 1;                     % Mean of initial human capital distribution
    nphi = 2;                     % Number of possible phi values (fixed productivity types)
    
    hlb = 0.0;                    % Human capital grid lower bound 
    hub = 7.0;                    % Human capital grid upper bound 
    nh = 3000;                    % Number of possible human capital values
    h_grid = linspace(hlb,hub,nh);% Create human capital grid
    
    simul_size = 20000;     % Number of individuals to follow in simulated sample
    nt = 400;               % Number of quarters to follow each individual - 100 yrs max - most will exit the model before this
    simul(simul_size,nt) = struct;
    % ======================= Create discretized h distribution ==================== %
    cdf_h = zeros(1,nh);          % CDF of initial h distribution
    pdf_h = zeros(1,nh);          % PDF of initial h distribution

    h_bin_size = (hub-hlb)/(nh-1);
    for ind_h = 1:nh
        if (ind_h ==1)
            bin_ub_now = h_grid(ind_h) + (h_bin_size)/2;
            cdf_h(ind_h) = normcdf(bin_ub_now,h_mu,h_sigma);
            pdf_h(ind_h) = normcdf(bin_ub_now,h_mu,h_sigma);
        else
            bin_ub_now = h_grid(ind_h) + (h_bin_size)/2;
            bin_ub_last = h_grid(ind_h-1) + (h_bin_size)/2;
            cdf_h(ind_h) = normcdf(bin_ub_now,h_mu,h_sigma);
            pdf_h(ind_h) = normcdf(bin_ub_now,h_mu,h_sigma) - normcdf(bin_ub_last,h_mu,h_sigma);
        end
    end
    
    npolicy=2;
    policy_urate = zeros(npolicy,1);
    policy_GDP = zeros(npolicy,1);
    policy_avg_hc = zeros(npolicy,1);
    policy_avg_wage = zeros(npolicy,1);
    policy_util = zeros(npolicy,1);
    policy_entrant = zeros(npolicy,1);
    
    for policy_ind =1:npolicy
        
        if policy_ind==2
            muA=1;
        end
    
        %% ============== Backward Loop over Ages to Get Value and Policy Functions ================ %%
        %---- Prespecify Value Functions
    
        U = zeros(nphi,nh);                % Value of unemployment in production phase                     
        U_hat  = zeros(nphi,nh);           % Value of unemployment in search phase - when searching for optimal job type
        U_hat_update  = zeros(nphi,nh);
    
        VN = zeros(nphi,nh);               % Value of employment in non-outsourced (N) job in production phase
        VO = zeros(nphi,nh);               % Value of employment in outsourced (O) job in production phase
        VN_hat = zeros(nphi,nh);           % Value of employment in N job in search phase
        VO_hat = zeros(nphi,nh);           % Value of employment in O job in search phase
        VN_hat_update = zeros(nphi,nh);
        VO_hat_update = zeros(nphi,nh);
        A_staffing = zeros(nphi,nh);
        %---- Prespecify Value Functions
        opt_ind_hp_N = ones(nphi,nh);     % Optimal index for next period h value (h') if in non-outsourced job
        opt_ind_hp_O = ones(nphi,nh);     % Opitmal index for next period h value (h') if in outsourced job
    
        GU = zeros(nphi,nh);              % Indicates optimal job type to search for if unemployed - 1=non-outsourced (N) and 2=outsourced (O)
        GN = zeros(nphi,nh);              % Indicates optimal job type to search for if employed in N job - 1=non-outsourced and 2=outsourced
        GO = zeros(nphi,nh);              % Indicates optimal job type to search for if employed in O job - 1=non-outsourced and 2=outsourced
    
        opt_p_U = zeros(nphi,nh);         % Optimal job-finding rate associated with optimal search choice when unemployed
        opt_p_N = zeros(nphi,nh);         % Optimal job-finding rate associated with optimal search choice when employed in N job
        opt_p_O = zeros(nphi,nh);         % Optimal job-finding rate associated with optimal search choice when employed in O job
    
        opt_x_U = zeros(nphi,nh);         % x contract value associated with optimal search choice when unemployed
        opt_x_N = zeros(nphi,nh);         % x contract value associated with optimal search choice when employed in N job
        opt_x_O = zeros(nphi,nh);         % x contract value associated with optimal search choice when employed in O job
    
        opt_eta_U = zeros(nphi,nh);       % eta (surplus rate) associated with optimal search choice when unemployed
        opt_eta_N = zeros(nphi,nh);       % eta (surplus rate) associated with optimal search choice when employed in N job 
        opt_eta_O = zeros(nphi,nh);       % eta (surplus rate) associated with optimal search choice when employed in O job
    
        opt_theta_U = zeros(nphi,nh);     % theta (market tightness) associated with optimal search choice when unemployed
        opt_theta_N = zeros(nphi,nh);     % theta (market tightness) associated with optimal search choice when employed in N job 
        opt_theta_O = zeros(nphi,nh);     % theta (market tightness) associated with optimal search choice when employed in O job
    
        % Find index for next period human capital (h') if the agent only lets human capital depreciate
        ind_hp_min = ones(1,nh);          % Minimum index for h' (index for h' if human capital only depreciates)
        for ind_h = 1:nh
            h = h_grid(ind_h);
            error_min = 100;
            for ind_hp = 1:nh
                hp = h_grid(ind_hp);       % hp is the value of h'
                error = abs(hp - h*(1-d)); % Find index that minimizes (hp-h*(1-d))
                if (error<error_min)
                    error_min = error;
                    ind_hp_min(ind_h) = ind_hp;
                end
            end
        end
    
        % Find index for next period human capital (h') if it is subject to the removal of some job-specific skill%
        ind_hp_js = ones(1,nh);              % Index for h' 
        for ind_h = 1:nh
            h = h_grid(ind_h);
            error_min = 100;
            for ind_hp = 1:nh
                hp = h_grid(ind_hp);         % hp is the value of h'
                error = abs(hp - h*(1-psi)); % Find index that minimizes (hp-h*(1-psi))
                if (error<error_min)
                    error_min = error;
                    ind_hp_js(ind_h) = ind_hp;
                end
            end
        end
    
        %--------------- Begin Value Function Iteration ----------------%
        vfi_error = 1; % Set vfi_error > vfi_toler
        while (vfi_error > vfi_toler)
    
            for ind_phi = 1:nphi          % Loop over phi types
                for ind_h = 1:nh          % Loop over current h values
                    h = h_grid(ind_h);                  % h = current human capital (h) value 
    
                    b_effect = b*z(ind_phi)*(h^alpha);  % Unemployment value
                    c_effect = c;                       % Search cost
    
                    % ---------- Solve hc Investment Decision ----------%
                    max_VN = -9999;
                    max_VO = -9999;
    
                    for ind_hp = ind_hp_min(ind_h):nh      % Loop over possible next period hp (h') values
                        hp = h_grid(ind_hp);
    
                        temp_VN = z(ind_phi)*(h^alpha) - ((hp - h*(1-d))*chi)^gamma ...
                            + beta*(1-delta)*VN_hat(ind_phi,ind_hp) + beta*delta*U_hat(ind_phi,ind_hp_js(ind_hp));
    
                        temp_VO = (1-muA)*(z(ind_phi)*(h^alpha)) - ((hp - h*(1-d))*chi)^gamma ...
                            + beta*(1-delta)*VO_hat(ind_phi,ind_hp) + beta*delta*U_hat(ind_phi,ind_hp_js(ind_hp));
    
                        if (temp_VN>max_VN)                % Save ind_hp where N job employment value is highest
                            max_VN = temp_VN;
                            VN(ind_phi,ind_h) = temp_VN;
                            opt_ind_hp_N(ind_phi,ind_h) = ind_hp;
                        end
                        if (temp_VO>max_VO)                % Save ind_hp where O job employment value is highest
                            max_VO = temp_VO;
                            VO(ind_phi,ind_h) = temp_VO;
                            opt_ind_hp_O(ind_phi,ind_h) = ind_hp;
                        end
    
                    end % end loop over hp options
                    U(ind_phi,ind_h) = b_effect + beta*U_hat(ind_phi,ind_hp_min(ind_h));
    
                    ind_hp = opt_ind_hp_O(ind_phi,ind_h);
                    A_staffing(ind_phi,ind_h) = muA*z(ind_phi)*(h^alpha) + (1-delta)*beta*(1-lambda*opt_p_O(ind_phi,ind_h))*A_staffing(ind_phi,ind_hp);
    
                    %---------------------------------------------------%
                    % ------- Solve Search Decision of Unemployed ----- %
                    xO = VO(ind_phi,ind_h) - c_effect;
                    etaO = (xO - U(ind_phi,ind_h))/(VO(ind_phi,ind_h) - U(ind_phi,ind_h));
                    t_Uhat_O = xO;
    
                    thetaN = (((VN(ind_phi,ind_h) - U(ind_phi,ind_h))/c_effect)^(ell/(1+ell)) - 1)^(1/ell);
                    if (~isreal(thetaN))
                        pN = 0;
                        xN = 0;
                        etaN = 0;
                    else
                        pN = thetaN/((1+thetaN^ell)^(1/ell));
                        qN = pN/thetaN;
                        xN = VN(ind_phi,ind_h) - (c_effect/qN);
                        etaN = (xN - U(ind_phi,ind_h))/(VN(ind_phi,ind_h) - U(ind_phi,ind_h));
                    end
                    t_Uhat_N = pN*xN + (1-pN)*U(ind_phi,ind_h);
    
                    % Now Consider which job type's value is highest
                    if (t_Uhat_N>t_Uhat_O)
                        GU(ind_phi,ind_h) = 1;
                        opt_p_U(ind_phi,ind_h) = pN;
                        U_hat_update(ind_phi,ind_h) = t_Uhat_N;
                        opt_x_U(ind_phi,ind_h) = xN;
                        opt_eta_U(ind_phi,ind_h) = etaN;
                        opt_theta_U(ind_phi,ind_h) = thetaN;
                    elseif (t_Uhat_O>t_Uhat_N)
                        GU(ind_phi,ind_h) = 2;
                        opt_p_U(ind_phi,ind_h) = 1;
                        U_hat_update(ind_phi,ind_h) = t_Uhat_O;
                        opt_x_U(ind_phi,ind_h) = xO;
                        opt_eta_U(ind_phi,ind_h) = etaO;
                        opt_theta_U(ind_phi,ind_h) = 1.0;
                    else
                        GU(ind_phi,ind_h) = 3;
                        opt_p_U(ind_phi,ind_h) = pN;
                        U_hat_update(ind_phi,ind_h) = t_Uhat_N;
                        opt_x_U(ind_phi,ind_h) = xN;
                        opt_eta_U(ind_phi,ind_h) = etaN;
                        opt_theta_U(ind_phi,ind_h) = thetaN;
                    end
    
                    %---------------------------------------------------%
                    % ------- Solve Search Decision of Employed --------%
                    xO = VO(ind_phi,ind_hp_js(ind_h)) - c_effect;
                    etaO = (xO - U(ind_phi,ind_hp_js(ind_h)))/(VO(ind_phi,ind_hp_js(ind_h)) - U(ind_phi,ind_hp_js(ind_h)));
    
                    t_VNhat_O = lambda*xO + (1-lambda)*VN(ind_phi,ind_h);
                    t_VOhat_O = lambda*xO + (1-lambda)*VO(ind_phi,ind_h);
    
                    thetaN_N = (((VN(ind_phi,ind_hp_js(ind_h)) - VN(ind_phi,ind_h))/c_effect)^(ell/(1+ell)) - 1)^(1/ell);
                    if (~isreal(thetaN_N))
                        pN_N = 0;
                        xN_N = 0;
                        etaN_N = 0;
                    else
                        pN_N = thetaN_N/((1+thetaN_N^ell)^(1/ell));
                        qN_N = pN_N/thetaN_N;
                        xN_N = VN(ind_phi,ind_hp_js(ind_h)) - (c_effect/qN_N);
                        etaN_N = (xN_N - U(ind_phi,ind_hp_js(ind_h)))/(VN(ind_phi,ind_hp_js(ind_h)) - U(ind_phi,ind_hp_js(ind_h)));
                    end
                    t_VNhat_N = lambda*pN_N*xN_N + (1-lambda*pN_N)*VN(ind_phi,ind_h);
    
                    thetaO_N = (((VN(ind_phi,ind_hp_js(ind_h)) - VO(ind_phi,ind_h))/c_effect)^(ell/(1+ell)) - 1)^(1/ell);
                    if (~isreal(thetaO_N))
                        pO_N = 0;
                        xO_N = 0;
                        etaO_N = 0;
                    else
                        pO_N = thetaO_N/((1+thetaO_N^ell)^(1/ell));
                        qO_N = pO_N/thetaO_N;
                        xO_N = VN(ind_phi,ind_hp_js(ind_h)) - (c_effect/qO_N);
                        etaO_N = (xO_N - U(ind_phi,ind_hp_js(ind_h)))/(VN(ind_phi,ind_hp_js(ind_h)) - U(ind_phi,ind_hp_js(ind_h)));
                    end  
                    t_VOhat_N = lambda*pO_N*xO_N + (1-lambda*pO_N)*VO(ind_phi,ind_h);
    
                    % For those Currently in N - Now Consider which job type is highest
                    if (t_VNhat_N>t_VNhat_O)
                        GN(ind_phi,ind_h) = 1;
                        opt_p_N(ind_phi,ind_h) = pN_N;
                        VN_hat_update(ind_phi,ind_h) = t_VNhat_N;
                        opt_x_N(ind_phi,ind_h) = xN_N;
                        opt_eta_N(ind_phi,ind_h) = etaN_N;
                        opt_theta_N(ind_phi,ind_h) = thetaN_N;
                    elseif (t_VNhat_O>t_VNhat_N)
                        GN(ind_phi,ind_h) = 2;
                        opt_p_N(ind_phi,ind_h) = 1;
                        VN_hat_update(ind_phi,ind_h) = t_VNhat_O;
                        opt_x_N(ind_phi,ind_h) = xO;
                        opt_eta_N(ind_phi,ind_h) = etaO;
                        opt_theta_N(ind_phi,ind_h) = 1.0;
                    else
                        GN(ind_phi,ind_h) = 3;
                        opt_p_N(ind_phi,ind_h) = pN_N;
                        VN_hat_update(ind_phi,ind_h) = t_VNhat_N;
                        opt_x_N(ind_phi,ind_h) = xN_N;
                        opt_eta_N(ind_phi,ind_h) = etaN_N;
                        opt_theta_N(ind_phi,ind_h) = thetaN_N;
                    end
    
                    % For those Currently in O - Now Consider which job type is highest
                    if (t_VOhat_N>t_VOhat_O)
                        GO(ind_phi,ind_h) = 1;
                        opt_p_O(ind_phi,ind_h) = pO_N;
                        VO_hat_update(ind_phi,ind_h) = t_VOhat_N;
                        opt_x_O(ind_phi,ind_h) = xO_N;
                        opt_eta_O(ind_phi,ind_h) = etaO_N;
                        opt_theta_O(ind_phi,ind_h) = thetaO_N;
                    elseif (t_VOhat_O>t_VOhat_N)
                        GO(ind_phi,ind_h) = 2;
                        opt_p_O(ind_phi,ind_h) = 1;
                        VO_hat_update(ind_phi,ind_h) = t_VOhat_O;
                        opt_x_O(ind_phi,ind_h) = xO;
                        opt_eta_O(ind_phi,ind_h) = etaO;
                        opt_theta_O(ind_phi,ind_h) = 1.0;
                    else 
                        GO(ind_phi,ind_h) = 3;
                        opt_p_O(ind_phi,ind_h) = pO_N;
                        VO_hat_update(ind_phi,ind_h) = t_VOhat_N;
                        opt_x_O(ind_phi,ind_h) = xO_N;
                        opt_eta_O(ind_phi,ind_h) = etaO_N;
                        opt_theta_O(ind_phi,ind_h) = thetaO_N;
                    end
    
                end % End loop over h values
            end % End loop over phi values
    
            error_U_hat = max(max(abs(U_hat_update(:,:) - U_hat(:,:))))/max(max(abs(U_hat_update(:,:))));
            error_VN_hat = max(max(abs(VN_hat_update(:,:) - VN_hat(:,:))))/max(max(abs(VN_hat_update(:,:))));
            error_VO_hat = max(max(abs(VO_hat_update(:,:) - VO_hat(:,:))))/max(max(abs(VO_hat_update(:,:))));
    
            vfi_error = max([error_U_hat,error_VN_hat,error_VO_hat]);
    
            U_hat(:,:) = U_hat_update(:,:);
            VN_hat(:,:) = VN_hat_update(:,:);
            VO_hat(:,:) = VO_hat_update(:,:);
    
        end % End while (vfi_error > vfi_toler)  
    
        %==================================================================================================================%
        %% ================================ Solve for Steady State Distribution ========================================= %%
        %==================================================================================================================%
    
        dist_U = zeros(nphi,nh);           % Steady state mass of unemployed agents
        dist_N = zeros(nphi,nh);           % Steady state mass of agents employed in non-outsourced job
        dist_O = zeros(nphi,nh);           % Steady state mass of agents employed in outsourced job
    
        dist_U_update = zeros(nphi,nh);
        dist_N_update = zeros(nphi,nh);
        dist_O_update = zeros(nphi,nh);
    
        dist_U(1,1) = pr_phi(1);           % Start with a guess regarding how the mass of agents is distributed
        dist_U(2,1) = pr_phi(2);
    
        dist_error = 1;
        while (dist_error>dist_toler)      % Iterate until distribution guess equals resulting distribution
    
            for ind_phi = 1:nphi           % Loop over phi types
                for ind_h = 1:nh           % Loop over current h values 
                    %=================================%
                    %========= New Entrants ==========%
                    %=================================%
                    % New Entrant - Finds Job
                    if (GU(ind_phi,ind_h)==1 || GU(ind_phi,ind_h)==3)     % Searches for N
                        dist_N_update(ind_phi,ind_h) = dist_N_update(ind_phi,ind_h)...
                            + opt_p_U(ind_phi,ind_h)*pr_phi(ind_phi)*nu*pdf_h(ind_h);
                    else                                                   % Searches for O
                        dist_O_update(ind_phi,ind_h) = dist_O_update(ind_phi,ind_h)...
                            + opt_p_U(ind_phi,ind_h)*pr_phi(ind_phi)*nu*pdf_h(ind_h);
                    end
                    % New Entrant - Does Not Find Job
                    dist_U_update(ind_phi,ind_h) = dist_U_update(ind_phi,ind_h) + (1-opt_p_U(ind_phi,ind_h))*pr_phi(ind_phi)*nu*pdf_h(ind_h);
    
                    %=================================%
                    %===== Previously Unemployed =====%
                    %=================================%
                    ind_hp = ind_hp_min(ind_h);             % New h value for the unemployed is what they get from allowing human capital to depreciate
    
                    % Unemployed - Finds a job
                    if (GU(ind_phi,ind_hp) == 1 || GU(ind_phi,ind_hp) == 3)     % Searches for N
                        dist_N_update(ind_phi,ind_hp) = dist_N_update(ind_phi,ind_hp)...
                            + opt_p_U(ind_phi,ind_hp)*dist_U(ind_phi,ind_h)*(1-nu);
                    else                                                         % Searches for O
                        dist_O_update(ind_phi,ind_hp) = dist_O_update(ind_phi,ind_hp)...
                            + opt_p_U(ind_phi,ind_hp)*dist_U(ind_phi,ind_h)*(1-nu);
                    end 
                    % Unemployed - Does not find a job
                    dist_U_update(ind_phi,ind_hp) = dist_U_update(ind_phi,ind_hp) + (1-opt_p_U(ind_phi,ind_hp))*dist_U(ind_phi,ind_h)*(1-nu);
    
                    %=================================%
                    %======= Previously in N =========%
                    %=================================%
                    ind_hp = opt_ind_hp_N(ind_phi,ind_h);   % New h value for employed in N job is the result of their optimal investment choice
    
                    % ----- N - No Separation Shock - Finds New Job
                    if (GN(ind_phi,ind_hp)==1 || GN(ind_phi,ind_hp)==3)      % N - Search OTJ for N
                        dist_N_update(ind_phi,ind_hp_js(ind_hp)) = dist_N_update(ind_phi,ind_hp_js(ind_hp))...
                            + lambda*opt_p_N(ind_phi,ind_hp)*(1-delta)*dist_N(ind_phi,ind_h)*(1-nu);
                    else                                                      % N - Search OTJ for O
                        dist_O_update(ind_phi,ind_hp_js(ind_hp)) = dist_O_update(ind_phi,ind_hp_js(ind_hp))...
                            + lambda*opt_p_N(ind_phi,ind_hp)*(1-delta)*dist_N(ind_phi,ind_h)*(1-nu);
                    end
    
                    % ----- N - No Separation Shock - Does Not Find New Job
                    dist_N_update(ind_phi,ind_hp) = dist_N_update(ind_phi,ind_hp)...
                        + (1-lambda*opt_p_N(ind_phi,ind_hp))*(1-delta)*dist_N(ind_phi,ind_h)*(1-nu);
    
                    % ----- N - Separation Shock - Finds Job
                    if (GU(ind_phi,ind_hp_js(ind_hp))==1 || GU(ind_phi,ind_hp_js(ind_hp))==3)   % N - Separation - Search for N
                        dist_N_update(ind_phi,ind_hp_js(ind_hp)) = dist_N_update(ind_phi,ind_hp_js(ind_hp))...
                            + opt_p_U(ind_phi,ind_hp_js(ind_hp))*delta*dist_N(ind_phi,ind_h)*(1-nu);
                    else                                                                        % N - Separation - Search for O
                        dist_O_update(ind_phi,ind_hp_js(ind_hp)) = dist_O_update(ind_phi,ind_hp_js(ind_hp))...
                            + opt_p_U(ind_phi,ind_hp_js(ind_hp))*delta*dist_N(ind_phi,ind_h)*(1-nu);
                    end
    
                    % ----- N - Separation Shock - Does Not Find New Job
                    dist_U_update(ind_phi,ind_hp_js(ind_hp)) = dist_U_update(ind_phi,ind_hp_js(ind_hp))...
                        + (1-opt_p_U(ind_phi,ind_hp_js(ind_hp)))*delta*dist_N(ind_phi,ind_h)*(1-nu);
    
                    %=================================%
                    %======= Previously in O =========%
                    %=================================%
                    ind_hp = opt_ind_hp_O(ind_phi,ind_h);   % New h value for employed in O job is the result of their optimal investment choice
    
                    % ----- O - No Separation Shock - Finds New Job
                    if (GO(ind_phi,ind_hp)==1 || GO(ind_phi,ind_hp)==3)      % O - Search OTJ for N
                        dist_N_update(ind_phi,ind_hp_js(ind_hp)) = dist_N_update(ind_phi,ind_hp_js(ind_hp))...
                            + lambda*opt_p_O(ind_phi,ind_hp)*(1-delta)*dist_O(ind_phi,ind_h)*(1-nu);
                    else                                                      % O - Search OTJ for O
                        dist_O_update(ind_phi,ind_hp_js(ind_hp)) = dist_O_update(ind_phi,ind_hp_js(ind_hp))...
                            + lambda*opt_p_O(ind_phi,ind_hp)*(1-delta)*dist_O(ind_phi,ind_h)*(1-nu);
                    end
                    % ----- O - No Separation Shock - Does Not Find New Job
                    dist_O_update(ind_phi,ind_hp) = dist_O_update(ind_phi,ind_hp)...
                        + (1-lambda*opt_p_O(ind_phi,ind_hp))*(1-delta)*dist_O(ind_phi,ind_h)*(1-nu);
    
                    % ----- O - Separation Shock - Finds Job
                    if (GU(ind_phi,ind_hp_js(ind_hp))==1 || GU(ind_phi,ind_hp_js(ind_hp))==3)  % O - Separation - Search for N
                        dist_N_update(ind_phi,ind_hp_js(ind_hp)) = dist_N_update(ind_phi,ind_hp_js(ind_hp))...
                            + opt_p_U(ind_phi,ind_hp_js(ind_hp))*delta*dist_O(ind_phi,ind_h)*(1-nu);
                    else                                                                       % O - Separation - Search for O
                        dist_O_update(ind_phi,ind_hp_js(ind_hp)) = dist_O_update(ind_phi,ind_hp_js(ind_hp))...
                            + opt_p_U(ind_phi,ind_hp_js(ind_hp))*delta*dist_O(ind_phi,ind_h)*(1-nu);
                    end
    
                    % ----- O - Separation Shock - Does Not Find New Job
                    dist_U_update(ind_phi,ind_hp_js(ind_hp)) = dist_U_update(ind_phi,ind_hp_js(ind_hp))...
                        + (1-opt_p_U(ind_phi,ind_hp_js(ind_hp)))*delta*dist_O(ind_phi,ind_h)*(1-nu);
    
                end % End loop over last period h value
            end % End loop over phi values
    
            error_dist_U = max(max(abs(dist_U_update(:,:) - dist_U(:,:))));
            error_dist_N = max(max(abs(dist_N_update(:,:) - dist_N(:,:))));
            error_dist_O = max(max(abs(dist_O_update(:,:) - dist_O(:,:))));
    
            dist_error = max([error_dist_U,error_dist_N,error_dist_O]);
    
            dist_U(:,:) = dist_U_update(:,:);
            dist_N(:,:) = dist_N_update(:,:);
            dist_O(:,:) = dist_O_update(:,:);
    
            dist_U_update(:,:) = 0;
            dist_N_update(:,:) = 0;
            dist_O_update(:,:) = 0;
    
    
        end % end while (dist_error>dist_toler)
    
        %% =============================================================================================================== %%
        %======================================== Simulation for Wage Moments and Results ==================================%
        %================================================================================================================= %%
    
        % Need to follow "respondents" from entrance because SS distribution does not save the eta values %
        rng(123456789);         % Set random number seed so results are exactly the same each time the code is ran   
    
        for id = 1:simul_size   % Loop over each individual in the sample

            for t = 1:nt        % Loop over each time period the individual is observed
    
                if (t==1)
                    % For new entrants, assign phi and h, then see if they are successful in their job search from unemployment %
                    % Assign starting h value
                    r1 = rand;
                    simul(id,t).h_ind = 1;  % If r1 is not bigger than any CDF values, assign to the lowest slot
                    for ind_h = 2:nh
                        if (r1>cdf_h(ind_h-1)) 
                            simul(id,t).h_ind = ind_h;
                        end 
                    end
                    ind_h = simul(id,t).h_ind;
                    h = h_grid(ind_h);
                    % Now assign phi type
                    r1 = rand;
                    if r1 < pr_phi(1)
                        simul(id,t).phi_ind = 1;
                    else
                        simul(id,t).phi_ind = 2;
                    end
                    ind_phi = simul(id,t).phi_ind;
                    % See if new entrant is successful in job search from unemployment - - If successful, save their wage %
                    r1 = rand;
                    if (r1 < opt_p_U(ind_phi,ind_h)) % Finds job
                        if (GU(ind_phi,ind_h)==1 || GU(ind_phi,ind_h)==3)        % Search for N
                            simul(id,t).stat = 1;                                % Status "stat" 1 = in non-outsourced job
                            simul(id,t).eta = opt_eta_U(ind_phi,ind_h);
                            eta = simul(id,t).eta;
    
                            ind_hp = opt_ind_hp_N(ind_phi,ind_h);                % This is the h' value they will choose for next period 
                            hp = h_grid(ind_hp);
                            simul(id,t).wage = eta*(VN(ind_phi,ind_h) - U(ind_phi,ind_hp_js(ind_h))) + U(ind_phi,ind_hp_js(ind_h)) + (chi*(hp - h*(1-d)))^gamma...
                                - beta*delta*U_hat(ind_phi,ind_hp_js(ind_hp)) - beta*(1-delta)*lambda*opt_p_N(ind_phi,ind_hp)*opt_x_N(ind_phi,ind_hp)...
                                - beta*(1-delta)*(1-lambda*opt_p_N(ind_phi,ind_hp))*(eta*VN(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp_js(ind_hp)));
    
                        else                                                      % Search for O
                            simul(id,t).stat = 2;                                 % Status "stat" 2 = in outsourced job 
                            simul(id,t).eta = opt_eta_U(ind_phi,ind_h);
                            eta = simul(id,t).eta;
    
                            ind_hp = opt_ind_hp_O(ind_phi,ind_h);                % This is the h' value they will choose for next period 
                            hp = h_grid(ind_hp);
                            simul(id,t).wage = eta*(VO(ind_phi,ind_h) - U(ind_phi,ind_hp_js(ind_h))) + U(ind_phi,ind_hp_js(ind_h)) + (chi*(hp - h*(1-d)))^gamma...
                                - beta*delta*U_hat(ind_phi,ind_hp_js(ind_hp)) - beta*(1-delta)*lambda*opt_p_O(ind_phi,ind_hp)*opt_x_O(ind_phi,ind_hp)...
                                - beta*(1-delta)*(1-lambda*opt_p_O(ind_phi,ind_hp))*(eta*VO(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp_js(ind_hp)));   
                        end   
                    else  % Does not find job
                        simul(id,t).stat = 0;                                    % Status "stat" 0 = unemployed  
                        simul(id,t).wage = 0.0;
                        simul(id,t).eta = 0.0;
                    end   % End if t=1 individual finds job or not
    
                else % so t>1
                    simul(id,t).phi_ind = simul(id,t-1).phi_ind;                 % Individual's phi value is permanent
                    ind_phi = simul(id,t).phi_ind;
    
                    if (simul(id,t-1).stat==999) % if dead                       % If the agent already exited the model they remain dead             
                        simul(id,t).stat=999;                                    % Status "stat" 999 = dead
                    else   % was not dead in last period
                        rnu = rand;
                        if (rnu<nu) % Dies at end of last period
                            simul(id,t).stat=999;
                        else        % Does NOT die at end of last period
    
                            if ((simul(id,t-1).stat)==0)        % If unemployed last period 
                                % Recall h from last period and apply any depreciation
                                ind_h = simul(id,t-1).h_ind;
                                simul(id,t).h_ind =  ind_hp_min(ind_h);
                                ind_h = simul(id,t).h_ind;
                                h = h_grid(ind_h);
    
                                % If unemployed last period - see if matches from unemployment !
                                r1 = rand;
                                if (r1 < opt_p_U(ind_phi,ind_h))    % Finds job
                                    if (GU(ind_phi,ind_h)==1 || GU(ind_phi,ind_h)==3)      % Search for N job
                                        simul(id,t).stat = 1; 
                                        simul(id,t).eta = opt_eta_U(ind_phi,ind_h);
                                        eta = simul(id,t).eta;
                                        ind_hp = opt_ind_hp_N(ind_phi,ind_h);
                                        hp = h_grid(ind_hp);
                                        simul(id,t).wage = eta*(VN(ind_phi,ind_h) - U(ind_phi,ind_hp_js(ind_h))) + U(ind_phi,ind_hp_js(ind_h)) + (chi*(hp - h*(1-d)))^gamma...
                                            - beta*delta*U_hat(ind_phi,ind_hp_js(ind_hp)) - beta*(1-delta)*lambda*opt_p_N(ind_phi,ind_hp)*opt_x_N(ind_phi,ind_hp)...
                                            - beta*(1-delta)*(1-lambda*opt_p_N(ind_phi,ind_hp))*(eta*VN(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp_js(ind_hp)));  
                                    else                                                    % Search for O job
                                        simul(id,t).stat = 2; 
                                        simul(id,t).eta = opt_eta_U(ind_phi,ind_h);
                                        eta = simul(id,t).eta;
                                        ind_hp = opt_ind_hp_O(ind_phi,ind_h);
                                        hp = h_grid(ind_hp);
                                        simul(id,t).wage = eta*(VO(ind_phi,ind_h) - U(ind_phi,ind_hp_js(ind_h))) + U(ind_phi,ind_hp_js(ind_h)) + (chi*(hp - h*(1-d)))^gamma...
                                            - beta*delta*U_hat(ind_phi,ind_hp_js(ind_hp)) - beta*(1-delta)*lambda*opt_p_O(ind_phi,ind_hp)*opt_x_O(ind_phi,ind_hp)...
                                            - beta*(1-delta)*(1-lambda*opt_p_O(ind_phi,ind_hp))*(eta*VO(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp_js(ind_hp)));  
                                    end
                                else                                 % Unemployed last period and does not find job%
                                    simul(id,t).stat = 0;
                                    simul(id,t).wage = 0.0;
                                    simul(id,t).eta = 0.0;
                                end
    
                            elseif ((simul(id,t-1).stat)==1)    % If in N last period  
                                % Recall h from last period and apply any investment
                                ind_h = simul(id,t-1).h_ind;
                                simul(id,t).h_ind = opt_ind_hp_N(ind_phi,ind_h);
                                ind_h = simul(id,t).h_ind;
                                h = h_grid(ind_h);
    
                                % If N last period - see if separation shock %
                                rdelta = rand;
                                if (rdelta < delta)  % N last period - separation shock    %
                                    % Apply job-specific skill loss
                                    simul(id,t).h_ind = ind_hp_js(ind_h);
                                    ind_h = simul(id,t).h_ind;
                                    h = h_grid(ind_h);
    
                                    % See if finds job from unemployment
                                    r1 = rand;
                                    if (r1 < opt_p_U(ind_phi,ind_h)) % N last period - separation shock - finds job from unemployment
                                        if (GU(ind_phi,ind_h)==1 || GU(ind_phi,ind_h)==3)      % Search for N
                                            simul(id,t).stat = 1; 
                                            simul(id,t).eta = opt_eta_U(ind_phi,ind_h);
                                            eta = simul(id,t).eta;
                                            ind_hp = opt_ind_hp_N(ind_phi,ind_h);
                                            hp = h_grid(ind_hp);
                                            simul(id,t).wage = eta*(VN(ind_phi,ind_h) - U(ind_phi,ind_hp_js(ind_h))) + U(ind_phi,ind_hp_js(ind_h)) + (chi*(hp - h*(1-d)))^gamma...
                                                - beta*delta*U_hat(ind_phi,ind_hp_js(ind_hp)) - beta*(1-delta)*lambda*opt_p_N(ind_phi,ind_hp)*opt_x_N(ind_phi,ind_hp)...
                                                - beta*(1-delta)*(1-lambda*opt_p_N(ind_phi,ind_hp))*(eta*VN(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp_js(ind_hp)));  
                                        else                                                    % Search for O
                                            simul(id,t).stat = 2; 
                                            simul(id,t).eta = opt_eta_U(ind_phi,ind_h);
                                            eta = simul(id,t).eta;
                                            ind_hp = opt_ind_hp_O(ind_phi,ind_h);
                                            hp = h_grid(ind_hp);
                                            simul(id,t).wage = eta*(VO(ind_phi,ind_h) - U(ind_phi,ind_hp_js(ind_h))) + U(ind_phi,ind_hp_js(ind_h)) + (chi*(hp - h*(1-d)))^gamma...
                                                - beta*delta*U_hat(ind_phi,ind_hp_js(ind_hp)) - beta*(1-delta)*lambda*opt_p_O(ind_phi,ind_hp)*opt_x_O(ind_phi,ind_hp)...
                                                - beta*(1-delta)*(1-lambda*opt_p_O(ind_phi,ind_hp))*(eta*VO(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp_js(ind_hp))); 
                                        end
                                    else                                 % N last period - separation shock - does find NOT job from unemployment
                                        simul(id,t).stat = 0;
                                        simul(id,t).wage = 0.0;
                                        simul(id,t).eta = 0.0;
                                    end
    
                                else                   % N last period - no separation shock %
                                    % See if successfully searchs OTJ
                                    r1 = rand;
                                    if (r1 < lambda*opt_p_N(ind_phi,ind_h)) % N last period - no separation shock - successfully searchs OTJ
    
    
                                        if (GN(ind_phi,ind_h)==1 || GN(ind_phi,ind_h)==3)             % N - Search OTJ for N
                                            simul(id,t).stat = 1; 
                                            simul(id,t).eta = opt_eta_N(ind_phi,ind_h);
                                            eta = simul(id,t).eta;
                                            % Apply job-specific skill loss
                                            simul(id,t).h_ind = ind_hp_js(ind_h);
                                            ind_h = simul(id,t).h_ind;
                                            h = h_grid(ind_h);    
    
                                            ind_hp = opt_ind_hp_N(ind_phi,ind_h);
                                            hp = h_grid(ind_hp);
                                            simul(id,t).wage = eta*(VN(ind_phi,ind_h) - U(ind_phi,ind_hp_js(ind_h))) + U(ind_phi,ind_hp_js(ind_h)) + (chi*(hp - h*(1-d)))^gamma...
                                                - beta*delta*U_hat(ind_phi,ind_hp_js(ind_hp)) - beta*(1-delta)*lambda*opt_p_N(ind_phi,ind_hp)*opt_x_N(ind_phi,ind_hp)...
                                                - beta*(1-delta)*(1-lambda*opt_p_N(ind_phi,ind_hp))*(eta*VN(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp_js(ind_hp))); 
                                        else                                                           % N - Search OTJ for O
                                            simul(id,t).stat = 2; 
                                            simul(id,t).eta = opt_eta_N(ind_phi,ind_h);
                                            eta = simul(id,t).eta;
                                            % Apply job-specific skill loss
                                            simul(id,t).h_ind = ind_hp_js(ind_h);
                                            ind_h = simul(id,t).h_ind;
                                            h = h_grid(ind_h); 
    
                                            ind_hp = opt_ind_hp_O(ind_phi,ind_h);
                                            hp = h_grid(ind_hp);
                                            simul(id,t).wage = eta*(VO(ind_phi,ind_h) - U(ind_phi,ind_hp_js(ind_h))) + U(ind_phi,ind_hp_js(ind_h)) + (chi*(hp - h*(1-d)))^gamma...
                                                - beta*delta*U_hat(ind_phi,ind_hp_js(ind_hp)) - beta*(1-delta)*lambda*opt_p_O(ind_phi,ind_hp)*opt_x_O(ind_phi,ind_hp)...
                                                - beta*(1-delta)*(1-lambda*opt_p_O(ind_phi,ind_hp))*(eta*VO(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp_js(ind_hp)));    
                                        end
                                    else                                         % N last period - no separation shock - does NOT successfully search OTJ 
                                        simul(id,t).stat = simul(id,t-1).stat;
                                        simul(id,t).eta = simul(id,t-1).eta;
                                        eta = simul(id,t).eta;
                                        ind_hp = opt_ind_hp_N(ind_phi,ind_h);
                                        hp = h_grid(ind_hp);
                                        simul(id,t).wage = eta*(VN(ind_phi,ind_h) - U(ind_phi,ind_hp_js(ind_h))) + U(ind_phi,ind_hp_js(ind_h)) + (chi*(hp - h*(1-d)))^gamma...
                                            - beta*delta*U_hat(ind_phi,ind_hp_js(ind_hp)) - beta*(1-delta)*lambda*opt_p_N(ind_phi,ind_hp)*opt_x_N(ind_phi,ind_hp)...
                                            - beta*(1-delta)*(1-lambda*opt_p_N(ind_phi,ind_hp))*(eta*VN(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp_js(ind_hp)));
                                    end
                                end % End if receives separation shock when N last period     
                            else                                                  % If in O last period  
                                % Recall h from last period and apply any investment
                                ind_h = simul(id,t-1).h_ind; 
                                simul(id,t).h_ind = opt_ind_hp_O(ind_phi,ind_h);
                                ind_h = simul(id,t).h_ind;
                                h = h_grid(ind_h);
    
                                % If O last period - see if separation shock %
                                rdelta = rand;
                                if (rdelta < delta) % O last period - separation shock    
                                    % Apply job-specific skill loss
                                    simul(id,t).h_ind = ind_hp_js(ind_h);
                                    ind_h = simul(id,t).h_ind;
                                    h = h_grid(ind_h); 
    
                                    % See if finds job from unemployment
                                    r1 = rand;
                                    if (r1 < opt_p_U(ind_phi,ind_h))    % N last period - separation shock - finds job from unemployment
                                        if (GU(ind_phi,ind_h)==1 || GU(ind_phi,ind_h)==3)       % Search for N
                                            simul(id,t).stat = 1; 
                                            simul(id,t).eta = opt_eta_U(ind_phi,ind_h);
                                            eta = simul(id,t).eta;
                                            ind_hp = opt_ind_hp_N(ind_phi,ind_h);
                                            hp = h_grid(ind_hp);
                                            simul(id,t).wage = eta*(VN(ind_phi,ind_h) - U(ind_phi,ind_hp_js(ind_h))) + U(ind_phi,ind_hp_js(ind_h)) + (chi*(hp - h*(1-d)))^gamma...
                                                - beta*delta*U_hat(ind_phi,ind_hp_js(ind_hp)) - beta*(1-delta)*lambda*opt_p_N(ind_phi,ind_hp)*opt_x_N(ind_phi,ind_hp)...
                                                - beta*(1-delta)*(1-lambda*opt_p_N(ind_phi,ind_hp))*(eta*VN(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp_js(ind_hp)));
                                        else                                                     % Search for O
                                            simul(id,t).stat = 2; 
                                            simul(id,t).eta = opt_eta_U(ind_phi,ind_h);
                                            eta = simul(id,t).eta;
                                            ind_hp = opt_ind_hp_O(ind_phi,ind_h);
                                            hp = h_grid(ind_hp);
                                            simul(id,t).wage = eta*(VO(ind_phi,ind_h) - U(ind_phi,ind_hp_js(ind_h))) + U(ind_phi,ind_hp_js(ind_h)) + (chi*(hp - h*(1-d)))^gamma...
                                                - beta*delta*U_hat(ind_phi,ind_hp_js(ind_hp)) - beta*(1-delta)*lambda*opt_p_O(ind_phi,ind_hp)*opt_x_O(ind_phi,ind_hp)...
                                                - beta*(1-delta)*(1-lambda*opt_p_O(ind_phi,ind_hp))*(eta*VO(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp_js(ind_hp))); 
                                        end
    
                                    else                                 % O last period - separation shock - does find NOT job from unemployment
                                        simul(id,t).stat = 0;
                                        simul(id,t).wage = 0.0;
                                        simul(id,t).eta = 0.0;
                                    end
                                else                   % O last period - no separation shock %
                                    %See if successfully searchs OTJ
                                    r1 = rand;
                                    if (r1 < lambda*opt_p_O(ind_phi,ind_h))   % O last period - no separation shock - successfully searchs OTJ
                                        if (GO(ind_phi,ind_h)==1 || GO(ind_phi,ind_h)==3)            % O - Search OTJ for N
                                            simul(id,t).stat = 1; 
                                            simul(id,t).eta = opt_eta_O(ind_phi,ind_h);
                                            eta = simul(id,t).eta;
                                            % Apply job-specific skill loss
                                            simul(id,t).h_ind = ind_hp_js(ind_h);
                                            ind_h = simul(id,t).h_ind;
                                            h = h_grid(ind_h); 
    
                                            ind_hp = opt_ind_hp_N(ind_phi,ind_h);
                                            hp = h_grid(ind_hp);
                                            simul(id,t).wage = eta*(VN(ind_phi,ind_h) - U(ind_phi,ind_hp_js(ind_h))) + U(ind_phi,ind_hp_js(ind_h)) + (chi*(hp - h*(1-d)))^gamma...
                                                - beta*delta*U_hat(ind_phi,ind_hp_js(ind_hp)) - beta*(1-delta)*lambda*opt_p_N(ind_phi,ind_hp)*opt_x_N(ind_phi,ind_hp)...
                                                - beta*(1-delta)*(1-lambda*opt_p_N(ind_phi,ind_hp))*(eta*VN(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp_js(ind_hp)));  
                                        else                                                          % O - Search OTJ for O
                                            simul(id,t).stat = 2; 
                                            simul(id,t).eta = opt_eta_O(ind_phi,ind_h);
                                            eta = simul(id,t).eta;
                                            % Apply job-specific skill loss
                                            simul(id,t).h_ind = ind_hp_js(ind_h);
                                            ind_h = simul(id,t).h_ind;
                                            h = h_grid(ind_h); 
    
                                            ind_hp = opt_ind_hp_O(ind_phi,ind_h);
                                            hp = h_grid(ind_hp);
                                            simul(id,t).wage = eta*(VO(ind_phi,ind_h) - U(ind_phi,ind_hp_js(ind_h))) + U(ind_phi,ind_hp_js(ind_h)) + (chi*(hp - h*(1-d)))^gamma...
                                                - beta*delta*U_hat(ind_phi,ind_hp_js(ind_hp)) - beta*(1-delta)*lambda*opt_p_O(ind_phi,ind_hp)*opt_x_O(ind_phi,ind_hp)...
                                                - beta*(1-delta)*(1-lambda*opt_p_O(ind_phi,ind_hp))*(eta*VO(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp_js(ind_hp)));  
                                        end
                                    else                                         % O last period - no separation shock - does NOT successfully searchs OTJ
                                        simul(id,t).stat = simul(id,t-1).stat;
                                        simul(id,t).eta = simul(id,t-1).eta;
                                        eta = simul(id,t).eta;
                                        ind_hp = opt_ind_hp_O(ind_phi,ind_h);  
                                        hp = h_grid(ind_hp);
                                        simul(id,t).wage = eta*(VO(ind_phi,ind_h) - U(ind_phi,ind_hp_js(ind_h))) + U(ind_phi,ind_hp_js(ind_h)) + (chi*(hp - h*(1-d)))^gamma...
                                            - beta*delta*U_hat(ind_phi,ind_hp_js(ind_hp)) - beta*(1-delta)*lambda*opt_p_O(ind_phi,ind_hp)*opt_x_O(ind_phi,ind_hp)...
                                            - beta*(1-delta)*(1-lambda*opt_p_O(ind_phi,ind_hp))*(eta*VO(ind_phi,ind_hp) + (1-eta)*U(ind_phi,ind_hp_js(ind_hp)));  
                                    end
                                end % End if receives separation shock when OL last period
                            end % End if for stat (when not dead or died before period start) last period
                        end         % End if dies at end of last period
                    end    % End if dead in last period
                end % End if t==1
            end % End loop over time periods in panel
        end % End loop over ids in simul_size
    
    
        count_O =0.0;
        count_N =0.0;
    
        count_E =0.0;
        count_U =0.0;
        for id=1:simul_size
            for t = 1:nt
                if (simul(id,t).stat==2)
                    count_O = count_O + 1.0;
                elseif (simul(id,t).stat==1)
                    count_N = count_N + 1.0;
                end
    
                if (simul(id,t).stat==0) 
                    count_U = count_U + 1.0;
                elseif (simul(id,t).stat==1 || simul(id,t).stat==2)
                    count_E = count_E + 1.0;
                end
            end
        end
    
        if policy_ind==1
            % ================================================================================================================ %
            %========================================== Compute Moments of Interest ===========================================%
            %================================================================================================================= %
    
            % ----------- m1: JTJ Probability
            jtj_rate = 0.0;
            m_sum = 0.0;
            for ind_phi=1:nphi
                for ind_h = 1:nh
    
                    ind_hp_N = opt_ind_hp_N(ind_phi,ind_h); 
                    ind_hp_O = opt_ind_hp_O(ind_phi,ind_h);
    
                    jtj_rate = jtj_rate + lambda*opt_p_N(ind_phi,ind_hp_N)*(1.0-delta)*dist_N(ind_phi,ind_h)*(1.0-nu);
                    jtj_rate = jtj_rate + lambda*opt_p_O(ind_phi,ind_hp_O)*(1.0-delta)*dist_O(ind_phi,ind_h)*(1.0-nu);
    
                    jtj_rate = jtj_rate + opt_p_U(ind_phi,ind_hp_js(ind_hp_N))*delta*dist_N(ind_phi,ind_h)*(1.0-nu);
                    jtj_rate = jtj_rate + opt_p_U(ind_phi,ind_hp_js(ind_hp_O))*delta*dist_O(ind_phi,ind_h)*(1.0-nu);
    
                    m_sum = m_sum + dist_N(ind_phi,ind_h)*(1.0-nu) + dist_O(ind_phi,ind_h)*(1.0-nu);
                end 
            end 
    
            if (m_sum>0.0) 
                jtj_rate = (jtj_rate/m_sum);
            else
                jtj_rate = 999.0;
            end
            m1 = jtj_rate;
    
            % ----------- m2: Unemployment Rate
            m2 = (sum(sum(dist_U(:,:)))/( sum(sum(dist_N(:,:))) + sum(sum(dist_O(:,:))) + sum(sum(dist_U(:,:))) ))*100.0;
    
            % ----------- m3: Total Percent Outsourced
            m3 = (( sum(sum(dist_O(:,:))) )/(sum(sum(dist_N(:,:))) + sum(sum(dist_O(:,:)))))*100.0;
    
            % ----------- m4: Ratio Outsourced With to Without Bachelor's Degree
            DO_phi1 = (( sum(dist_O(1,:)) )/( sum(dist_O(1,:)) + sum(dist_N(1,:)) ))*100.0;
            DO_phi2 = (( sum(dist_O(2,:)) )/( sum(dist_O(2,:)) + sum(dist_N(2,:)) ))*100.0;
            m4 = DO_phi1/DO_phi2;
    
            % ----------- m5: Avg. Quarterly Separation Rate
            delta_rate = 0.0;
            m_sum = 0.0;
    
            for ind_phi=1:nphi
                for ind_h = 1:nh
    
                    ind_hp_N = opt_ind_hp_N(ind_phi,ind_h); 
                    ind_hp_O = opt_ind_hp_O(ind_phi,ind_h);    
    
                    delta_rate = delta_rate + (1.0-opt_p_U(ind_phi,ind_hp_js(ind_hp_N)))*delta*dist_N(ind_phi,ind_h)*(1.0-nu);
                    delta_rate = delta_rate + (1.0-opt_p_U(ind_phi,ind_hp_js(ind_hp_O)))*delta*dist_O(ind_phi,ind_h)*(1.0-nu);
    
                    m_sum = m_sum + dist_N(ind_phi,ind_h)*(1.0-nu) + dist_O(ind_phi,ind_h)*(1.0-nu);
                end 
            end 
            m5 = delta_rate/m_sum;
    
            % ----------- m6: Average Wage Loss After Unemployment Spell   
            wg_uspell = 0;
            m_sum = 0;
            for id = 1:simul_size
                stat_u = 0;
                for t = 1:nt-2 % Need to observe employed wage, unemployment, and then employed wage again
    
                    if (simul(id,t).stat~=999)
                        if (stat_u==1 && simul(id,t).stat~=0) % Was Unemployed but found Employment
                            w_after = simul(id,t).wage;
                            wg_uspell = wg_uspell + (log(w_after) - log(w_before))*100;
                            m_sum = m_sum +1;
                            stat_u = 0;
                        end
    
                        if (simul(id,t).stat~=0 && simul(id,t+1).stat==0) % Employed in t and unemployed in t+1
                            w_before = simul(id,t).wage;
           
                            stat_u = 1;
                        end 
                    end % End if has not exited the model (if stat~=999)
                end
            end
            wg_uspell = wg_uspell/m_sum;
            m6 = wg_uspell;
    
            % ----------- m7: Wage Dispersion p90/p50
            count_wages = 0;
            for id = 1:simul_size
                for t = 1:nt       
                    if (simul(id,t).stat~=0 && simul(id,t).stat~=999)
                        count_wages = count_wages + 1;
                    end    
                end
            end
    
            % --- First save sample of wages
            wage_sim = zeros(1,count_wages);
            count_wages = 0;
            for id = 1:simul_size
                for t = 1:nt       
                    if (simul(id,t).stat~=0 && simul(id,t).stat~=999)
                        count_wages = count_wages + 1;
                        wage_sim(count_wages) = simul(id,t).wage;
                    end    
                end
            end
    
    
            p90_wage = prctile(wage_sim, 90);
            p50_wage = prctile(wage_sim, 50);
    
            m7 = p90_wage/p50_wage;
    
            %----------- m8: Average Annual Real Wage Growth
            % average wage growth over 1 year
            wg_sum1 = 0.0;
            m_sum = 0.0;
            for id = 1:simul_size
                for t = 1:(nt-4)
                    if (simul(id,t).stat~=0 && simul(id,t+4).stat~=0 && simul(id,t).stat~=999 && simul(id,t+4).stat~=999) 
                        w1 = simul(id,t).wage;
                        w5 = simul(id,t+4).wage; 
                        wg_sum1 = wg_sum1 + ((w5-w1)/w1)*100.0;
                        m_sum = m_sum + 1.0;
                    end 
                end 
            end 
            wg_sum1 = wg_sum1/m_sum;
    
            m8 = wg_sum1;
    
            % ----------- m9: College Wage Premium
            avg_wage_phi1 = 0.0;
            avg_wage_phi2 = 0.0;
            sum_phi1= 0.0;
            sum_phi2= 0.0;
    
            for id = 1:simul_size
                for t = 1:nt
    
                    if (simul(id,t).stat~=0 && simul(id,t).stat~=999)   % If Not Unemployed or dead
                        if (simul(id,t).phi_ind == 1) 
                            avg_wage_phi1 = avg_wage_phi1 + simul(id,t).wage;
                            sum_phi1 = sum_phi1 + 1.0;
                        else
                            avg_wage_phi2 = avg_wage_phi2 + simul(id,t).wage;
                            sum_phi2 = sum_phi2 + 1.0;
                        end 
                    end 
                end 
            end 
    
            avg_wage_phi1 = avg_wage_phi1/sum_phi1;
            avg_wage_phi2 = avg_wage_phi2/sum_phi2;
    
            m9 = avg_wage_phi2/avg_wage_phi1;
    
            %---------- m10 - Unemployment Rate Among Age 20
            count_E =0.0;
            count_U =0.0;
            for id = 1:simul_size
                for t = 1:1 % age = 20
                    if (simul(id,t).stat~=999) 
    
                        if (simul(id,t).stat==0)  
                            count_U = count_U + 1.0;
                        elseif (simul(id,t).stat==1 || simul(id,t).stat==2) 
                            count_E = count_E + 1.0;
                        end 
                    end 
                end 
            end
            m10 = (count_U/(count_E + count_U))*100.0;
    
    
    
            mm = [m1;m2;m3;m4;m5;m6;m7;m8;m9;m10];

            %--------------- Figures: Human Capital Investment and Distributions -----------------%
            opt_hp_N = h_grid(opt_ind_hp_N);
    
            
            % HC distributions
            cdf_N= zeros(2,nh);
            cdf_O= zeros(2,nh);
            for ind_h =1:nh
                for ind_phi=1:nphi
                    if (ind_h==1)
                        cdf_N(ind_phi,ind_h) = (dist_N(ind_phi,ind_h))/( sum(dist_N(ind_phi,:)));
                        cdf_O(ind_phi,ind_h) = (dist_O(ind_phi,ind_h))/( sum(dist_O(ind_phi,:)));
                    else
                        cdf_N(ind_phi,ind_h) = cdf_N(ind_phi,ind_h-1) + (dist_N(ind_phi,ind_h))/( sum(dist_N(ind_phi,:)) );
                        cdf_O(ind_phi,ind_h) = cdf_O(ind_phi,ind_h-1) + (dist_O(ind_phi,ind_h))/( sum(dist_O(ind_phi,:)) );
                    end
                end
            end
    
            %----------------  Optimal Investment Rule: h'-h ----------------%
            opt_hp_O = h_grid(opt_ind_hp_O);
            h_diff_N(1,:) = opt_hp_N(1,:) - h_grid(:)';
            h_diff_N(2,:) = opt_hp_N(2,:) - h_grid(:)';
    
            h_diff_O(1,:) = opt_hp_O(1,:) - h_grid(:)';
            h_diff_O(2,:) = opt_hp_O(2,:) - h_grid(:)';
    
            % N vs O - College and Non-College
            avg_diff_N = zeros(2,12);
            avg_diff_O = zeros(2,12);
            for ind_phi = 1:nphi
                avg_diff_N(ind_phi,1) = mean(h_diff_N(ind_phi,1:108));        %0-0.25
                avg_diff_N(ind_phi,2) = mean(h_diff_N(ind_phi,109:215));      %0.25-0.5
                avg_diff_N(ind_phi,3) = mean(h_diff_N(ind_phi,216:322));      %0.5-0.75
                avg_diff_N(ind_phi,4) = mean(h_diff_N(ind_phi,323:429));      %0.75-1
    
                avg_diff_N(ind_phi,5) = mean(h_diff_N(ind_phi,430:536));      %1-1.25
                avg_diff_N(ind_phi,6) = mean(h_diff_N(ind_phi,537:643));      %1.25-1.5
                avg_diff_N(ind_phi,7) = mean(h_diff_N(ind_phi,644:750));      %1.5-1.75
                avg_diff_N(ind_phi,8) = mean(h_diff_N(ind_phi,751:857));      %1.75-2
    
                avg_diff_N(ind_phi,9) = mean(h_diff_N(ind_phi,858:964));      %2-2.25
                avg_diff_N(ind_phi,10) = mean(h_diff_N(ind_phi,965:1072));    %2.25-2.5
                avg_diff_N(ind_phi,11) = mean(h_diff_N(ind_phi,1073:1179));   %2.5-2.75
                avg_diff_N(ind_phi,12) = mean(h_diff_N(ind_phi,1180:1286));   %2.75-3
    
                avg_diff_O(ind_phi,1) = mean(h_diff_O(ind_phi,1:108));        %0-0.25
                avg_diff_O(ind_phi,2) = mean(h_diff_O(ind_phi,109:215));      %0.25-0.5
                avg_diff_O(ind_phi,3) = mean(h_diff_O(ind_phi,216:322));      %0.5-0.75
                avg_diff_O(ind_phi,4) = mean(h_diff_O(ind_phi,323:429));      %0.75-1
    
                avg_diff_O(ind_phi,5) = mean(h_diff_O(ind_phi,430:536));      %1-1.25
                avg_diff_O(ind_phi,6) = mean(h_diff_O(ind_phi,537:643));      %1.25-1.5
                avg_diff_O(ind_phi,7) = mean(h_diff_O(ind_phi,644:750));      %1.5-1.75
                avg_diff_O(ind_phi,8) = mean(h_diff_O(ind_phi,751:857));      %1.75-2
    
                avg_diff_O(ind_phi,9) = mean(h_diff_O(ind_phi,858:964));      %2-2.25
                avg_diff_O(ind_phi,10) = mean(h_diff_O(ind_phi,965:1072));    %2.25-2.5
                avg_diff_O(ind_phi,11) = mean(h_diff_O(ind_phi,1073:1179));   %2.5-2.75
                avg_diff_O(ind_phi,12) = mean(h_diff_O(ind_phi,1180:1286));   %2.75-3
    
            end
    
    
        end % end if policy_ind==1
        %=============== Estimate Aggregate Effects of Removing Outsourcing Option =================%
        % Compute and save moments of interest for comparison
        %-- Unemployment rate
        policy_urate(policy_ind) = (sum(sum(dist_U(:,:)))/( sum(sum(dist_N(:,:))) + sum(sum(dist_O(:,:))) + sum(sum(dist_U(:,:))) ))*100;
    
        %----------- Output from Employment (GDP)
        policy_GDP(policy_ind) = 0;
        for ind_phi = 1:nphi
            for ind_h = 1:nh
                h = h_grid(ind_h);
                policy_GDP(policy_ind) = policy_GDP(policy_ind) + (z(ind_phi)*h^alpha)*(dist_N(ind_phi,ind_h) + dist_O(ind_phi,ind_h));
            end
        end
    
        %---------- Average Human Capital
        h_sum_total = 0;
        policy_avg_hc(policy_ind) = 0;
        for ind_phi = 1:nphi
            for ind_h = 1:nh
                h = h_grid(ind_h);
    
                policy_avg_hc(policy_ind) = policy_avg_hc(policy_ind) + h*(dist_U(ind_phi,ind_h) + dist_N(ind_phi,ind_h) + dist_O(ind_phi,ind_h));
                h_sum_total = h_sum_total + dist_U(ind_phi,ind_h) + dist_N(ind_phi,ind_h) + dist_O(ind_phi,ind_h);
            end
        end
        policy_avg_hc(policy_ind) = policy_avg_hc(policy_ind)/h_sum_total;
    
    
        %-------- Average Wage
        count_wage_total = 0;
        policy_avg_wage(policy_ind) = 0;
        for id = 1:simul_size
            for t = 1:nt
                if (simul(id,t).stat~=999 && simul(id,t).stat~=0) % If employed
                    policy_avg_wage(policy_ind) = policy_avg_wage(policy_ind) + simul(id,t).wage;
                    count_wage_total = count_wage_total + 1;
    
                end
            end
        end
        policy_avg_wage(policy_ind) = policy_avg_wage(policy_ind)/count_wage_total;
    
        %--------- Total Utility
        util = 0;
        for ind_phi = 1:nphi
            for ind_h = 1:nh
                h = h_grid(ind_h);
                hp_N = h_grid(opt_ind_hp_N(ind_phi,ind_h));
                hp_O = h_grid(opt_ind_hp_O(ind_phi,ind_h));
                % Compute production minus investment cost
                util = util + dist_U(ind_phi,ind_h)*(b*z(ind_phi)*(h^alpha))...
                    + z(ind_phi)*(h^alpha)*(dist_N(ind_phi,ind_h) + dist_O(ind_phi,ind_h))...
                    - (((hp_N - h*(1-d))*chi)^gamma)*dist_N(ind_phi,ind_h) - (((hp_O - h*(1-d))*chi)^gamma)*dist_O(ind_phi,ind_h);
                % Now subtract out vacancy posting costs
                if (GU(ind_phi,ind_h)==1 || GU(ind_phi,ind_h)==3) % Search in unemployment for non-outsourced job
                    if (~isnan(opt_theta_U(ind_phi,ind_h)) && isreal(opt_theta_U(ind_phi,ind_h)))
                        util = util - dist_U(ind_phi,ind_h)*opt_theta_U(ind_phi,ind_h)*c;
                    end
                else                                              % Search in unemployment for outsourced job
                    if (~isnan(opt_theta_U(ind_phi,ind_h)) && isreal(opt_theta_U(ind_phi,ind_h)))
                        util = util - dist_U(ind_phi,ind_h)*opt_theta_U(ind_phi,ind_h)*(c + A_staffing(ind_phi,ind_h));
                    end
                end
    
                if (GN(ind_phi,ind_h)==1 || GN(ind_phi,ind_h)==3) % Search from non-outsourced job for non-outsourced job
                    if (~isnan(opt_theta_N(ind_phi,ind_h)) && isreal(opt_theta_N(ind_phi,ind_h)))
                        util = util - dist_N(ind_phi,ind_h)*opt_theta_N(ind_phi,ind_h)*lambda*c;
                    end
                else                                              % Search from outsourced job for non-outsourced job
                    if (~isnan(opt_theta_N(ind_phi,ind_h))&& isreal(opt_theta_N(ind_phi,ind_h)))
                        util = util - dist_N(ind_phi,ind_h)*opt_theta_N(ind_phi,ind_h)*lambda*(c + A_staffing(ind_phi,ind_h));
                    end
                end
    
                if (GO(ind_phi,ind_h)==1 || GO(ind_phi,ind_h)==3) % Search from outsourced job for non-outsourced job
                    if (~isnan(opt_theta_O(ind_phi,ind_h)) && isreal(opt_theta_O(ind_phi,ind_h)))
                        util = util - dist_O(ind_phi,ind_h)*opt_theta_O(ind_phi,ind_h)*lambda*c;
                    end
                else
                    if (~isnan(opt_theta_O(ind_phi,ind_h)) && isreal(opt_theta_O(ind_phi,ind_h)))
                        util = util - dist_O(ind_phi,ind_h)*opt_theta_O(ind_phi,ind_h)*lambda*(c + A_staffing(ind_phi,ind_h));
                    end
                end
            end
        end
        policy_util(policy_ind) = util;
    
        %---------- New Entrant Value
        policy_entrant(policy_ind) = pr_phi(1)*U_hat(1,1) + pr_phi(2)*U_hat(2,1);
    
    
    end
    
    ppc_urate = policy_urate(2)-policy_urate(1);
    pc_GDP = ((policy_GDP(2) - policy_GDP(1))/policy_GDP(1))*100;
    pc_avg_hc = ((policy_avg_hc(2) - policy_avg_hc(1))/policy_avg_hc(1))*100;
    pc_avg_wage = ((policy_avg_wage(2) - policy_avg_wage(1))/policy_avg_wage(1))*100;
    pc_total_util = ((policy_util(2) - policy_util(1))/policy_util(1))*100;
    pc_entrant = ((policy_entrant(2) - policy_entrant(1))/policy_entrant(1))*100;
    
    policy_change_table = [ppc_urate; pc_GDP; pc_avg_hc; pc_avg_wage; pc_total_util; pc_entrant];

end % end function
    
    